Generation of third-harmonic spin oscillation from strong spin precession induced by terahertz magnetic near fields

The ability to drive a spin system to state far from the equilibrium is indispensable for investigating spin structures of antiferromagnets and their functional nonlinearities for spintronics. While optical methods have been considered for spin excitation, terahertz (THz) pulses appear to be a more convenient means of direct spin excitation without requiring coupling between spins and orbitals or phonons. However, room-temperature responses are usually limited to small deviations from the equilibrium state because of the relatively weak THz magnetic fields in common approaches. Here, we studied the magnetization dynamics in a HoFeO3 crystal at room temperature. A custom-made spiral-shaped microstructure was used to locally generate a strong multicycle THz magnetic near field perpendicular to the crystal surface; the maximum magnetic field amplitude of about 2 T was achieved. The observed time-resolved change in the Faraday ellipticity clearly showed second- and third-order harmonics of the magnetization oscillation and an asymmetric oscillation behaviour. Not only the ferromagnetic vector M but also the antiferromagnetic vector L plays an important role in the nonlinear dynamics of spin systems far from equilibrium.

Since their first use a decade ago, THz pulses have been considered a more efficient means of directly exciting spin waves (magnons) in antiferromagnets for spin control [17][18][19][20][21][22][23][24][25] . One of the distinct signatures of nonlinear spin dynamics is the generation of second-harmonic (SH) or third-harmonic (TH) oscillations. However, because the maximum peak amplitude of a single-cycle THz pulse usually reaches only about 0.1 T and its spectral density at the spin resonance frequency is also limited, the harmonic spin oscillation is limited to the second order in common approaches 26 . Although it has been shown that a THz electric field enhanced by an antenna structure allows us to rotate the magnetization at low temperatures 27 , the resulting spin dynamics showed no higher-order harmonic oscillation near the critical temperature of spin reorientation. Furthermore, the spin dynamics of antiferromagnets are described by two sublattices, and therefore both the ferromagnetic vector M and the antiferromagnetic vector L may contribute to the nonlinear responses 8 . However, the observation of a sizable L at room temperature requires strong excitation. Even a THz magnetic field enhanced by a split ring resonator was not able to reveal obvious nonlinear properties besides a redshift in the magnon frequency; 28 in other words, the excitation was still too weak. Such experimental results can be sufficiently described by the dynamics of M without the need of considering L. Thus, it has remained unclear how M and L contribute to the generation of higher-order harmonic spin oscillations in the strong excitation regime and nonlinear magnetization changes in general. In this study, we observed a large spin-precession amplitude that generates the SH and TH of spin motion in the canted antiferromagnet HoFeO 3 by using a large THz magnetic near field, and clarified their relation with the dynamics of the ferromagnetic vector M but also the antiferromagnetic vector L.

Results
Magnetic field enhancement in a spiral-shaped microstructure The measurement geometry is presented in Fig. 1a (Methods and Supplementary Section I). The sample (shown in blue) was a 52-μm-thick c-cut HoFeO 3 crystal with an intrinsic quasi-antiferromagnetic (q-AF) mode at ν AF = 0.58 THz. The linearly polarized electric field of a THz pulse (red thick arrow) was coupled to the long triangular tail of a metallic spiral-shaped microstructure on the surface of the sample (shown in yellow). The original waveform of the THz pulse is shown by the gray curve in Fig. 1b 29 . To remove the field components that are not needed to excite the q-AF mode, a low-pass filter was inserted in front of the sample (cut-off frequency: 0.68 THz; the black curve in Fig. 1b is the obtained time trace). An image of the fabricated structure is shown in Fig. 1c. The THz pulse-induced current flows into the spiral of the microstructure, which enhances the THz magnetic field. The magnetic field at the center of spiral structure is about 200 times larger than that of the incident THz magnetic field at the resonance frequency ν c of the structure [=0.54 THz, which was determined by a finite-difference time-domain (FDTD) simulation]. Compared with the case of using a split ring resonator 28 , a three times stronger and spatially smoother magnetic field can be realized by using this structure, as shown at the bottom of Fig. 1c (Supplementary Sections II-IV). Accordingly, we were able to generate magnetic field amplitudes up to 2.1 T even with the filtered THz pulse. The strong magnetic field along the z-axis (see the black curve in Fig. 1d) exerts a Zeeman torque on the spins in HoFeO 3 and thus induces a change in the net magnetization along the zaxis, ΔM z .

Asymmetric changes in the Faraday ellipticity
As shown by the thin red arrows in Fig. 1a, a linearly polarized nearinfrared (NIR) probe pulse for the Faraday ellipticity measurements is incident from the back of the HoFeO 3 crystal and focused at the center of the spiral. The measurement is based on the fact that a change in the magnetization leads to a change in the polarization ellipticity angle of the probe light that passes through the crystal (Supplementary Section I). The three temporal profiles of the change in the ellipticity angle (Δη) in Fig. 1d for THz electric field amplitudes equal to E THz = 0.8, 0.55, and 0.2 MV cm −1 correspond to the data for B z = 2.1, 1.4, and 0.5 T at the sample surface, respectively. We can see      that the curve for E THz = 0.8 MV cm −1 exhibits a distorted sinusoidal oscillation with an asymmetric amplitude distribution; i.e., the amplitudes of the signals with Δη > 0 are larger than those of the signals with Δη < 0.
The asymmetric waveform of Δη for 0.8 MV cm −1 observed over the course of about 30 ps is a fingerprint of a large magnetization change induced by the high magnetic field. This data is not a result of the THz Kerr effect, because the Kerr effect is attributed to the diagonal elements of the third-order nonlinear susceptibility χ ð3Þ and thus should not depend on the sign of the magnetization, but our measured Δη depends on the sign as shown in Fig. 2 and Supplementary Section V. We can also exclude thermally induced changes in the spin structure, such as demagnetization or spin reorientation, because their recovery times (to reach the equilibrium state) are much longer than tens of picoseconds 30 .

Data analysis using ferromagnetic and antiferromagnetic vectors
To explain the waveform of Δη and to evaluate the magnetization change achieved by our method, we numerically solved the Landau-Lifshitz-Gilbert (LLG) equation and the propagation equation of the probe pulse (Supplementary Sections VI and VII). The former equation was used to obtain the dynamics of the sublattice magnetizations m i with i = 1, 2, and the latter equation was used to calculate the Faraday ellipticity change occurring in the sample. The propagation equation includes the following magnetization-dependent permittivity tensor for HoFeO 3 : 31 Here, the parameters ε 0 and σ are magnetization-independent terms, whereas the off-diagonal term κ is a magnetization-dependent term. σ describes the strength of the birefringence. Regarding κ, in our calculation, the contributions from both the ferromagnetic vector M = m 1 + m 2 and the antiferromagnetic vector L = m 1 − m 2 are taken into account (note that L can only be ignored in the case of small deflection angles): 31 where ΔM z (ΔL x ) is the temporal variation of the z-axis (x-axis) component of M (L). M 0 and L 0 are the initial amplitudes of these vectors. f and g represent phenomenological constants. The calculated waveforms for different magnetic field strengths are shown by the open circles in Fig. 1d, which reveals that our model reproduces the experimental data (solid curves) well. Our calculation also reproduces the absolute values of η for different signs of the initial magnetization, ±M 0 (Fig. 2).

Interpretation of the asymmetric temporal profiles
The good agreement between the experiment and the calculation allows us to understand the origin of the asymmetric temporal profiles of Δη: they originate from spin dynamics with a large deflection angle. The calculated dynamics of ΔM z /M 0 and ΔL x /L 0 are presented in Fig. 3a, b, respectively. In stark contrast to the behavior of ΔM z /M 0 , ΔL x /L 0 shows a vertically asymmetrical behavior. This difference can be graphically explained by the motion of L: The red (blue) trajectories in Fig. 3c show a condition where the spins are strongly (weakly) excited and deviate from the x-axis by a large (small) maximum deflection angle θ. While M oscillates and its amplitude changes, L rotates around the zaxis by ±θ (with an almost constant | L| due to the small canting angles of m i ; see Supplementary Section VIII). The rotation of L is shown in Fig. 3d (the top view of Fig. 3c). As shown by the calculated results in Fig. 3b, the x component of L always decreases as the deflection angle increases and ΔL x oscillates with twice the frequency of ΔM z . The difference in the oscillation frequency is explained in Fig. 3d: the change in the x component of L always proceeds along the same path (x = 2 → 2-| ΔL x |→ 2), independently of the path of the arrowhead of the vector L projected on the y-axis (either y = 0 → +|ΔL y |→ 0 or 0 → -|ΔL y |→ 0). The field-strength dependences of their maximum values (ΔM max z =M 0 , ΔL max x =L 0 , and θ) are presented in Fig. 3e, where θ almost linearly depends on E THz and ΔM max z =M 0 increases approximately linearly with θ. However, because ΔL max x =L 0 exhibits a cosine relationship with θ, it is much less sensitive to changes in θ as long as θ is small. Thus, only if the spin precession is driven far away from the equilibrium state, a considerable change in ΔL x /L 0 appears and contributes to Δη, causing the asymmetry shown in Fig. 1d. The data points indicated by the black arrows in Fig. 3e represent the values calculated for the electric fields used in the experiment. From the results, we determine θ ≈ 40°, ΔM max z =M 0 ≈ 0.9 and ΔL max x =L 0 ≈ 0.2 at the surface of HoFeO 3 for E THz = 0.8 MV cm −1 . In the experiment in ref. 28

Higher-order harmonic spectra in the time-frequency domain
In addition to the observed asymmetric change in magnetization, the generation of an SH or TH spin oscillation is a hallmark of a large magnetization amplitude far away from the equilibrium value M 0 . As shown in Fig. 4a, the Fourier transform spectra of the experimentally obtained Δη(t) reveal the nonlinearity of the spin motion caused by the enhanced THz magnetic field. In the spectrum for E THz = 0.2 MV cm −1 (blue spectrum), there is a fundamental peak at ν = 0.58 THz, which is equal to the q-AF mode frequency ν AF , and another peak lies at 1.16 THz, which precisely matches the SH of ν AF . In the spectrum for E THz = 0.8 MV cm −1 (red shaded area in Fig. 4a), there is a peak corresponding to the TH, in addition to the fundamental and SH peaks. However, compared with the spectrum measured at weaker fields (blue shaded area in Fig. 4a), the center frequencies of the SH (1.06 THz) and TH (1.58 THz) peaks are slightly lower than 2ν AF (1.16 THz) and 3ν AF (1.74 THz), and the three peaks exhibit an obvious broadening towards the low-frequency side.
To elaborate on the observed behavior of the higher-order harmonic oscillations, we examine the time-resolved Fourier transform spectra of Δη(t) for E THz = 0.8 and 0.2 MV cm −1 in Fig. 4b. While the maximum amplitude of the fundamental peak becomes stronger as E THz increases, the rise time and decay time become shorter. To focus on the signal duration, the temporal dynamics of the normalized fundamental, SH, and TH signals obtained by integrating the peak area are shown in Fig. 4c.(the data without normalization are shown in Supplementary Section IX). The signals of the SH and TH peaks reach their maxima at almost the same time as the fundamental peak, but they decay faster, verifying that the SH and TH are indeed generated only when the fundamental magnetization change is large. The spectral position of the fundamental peak as a function of time is shown in Fig. 4d.at early times before 20 ps, the oscillation peak frequencies notably deviate from ν AF and the oscillation peak frequency even reaches ν red = 0.53 THz for E THz = 0.8 MV cm −1 . This redshift is a result of the intrinsic nonlinear properties of the antiferromagnet in addition to the effect of the forced oscillation of the spin in the THz magnetic field at ν c = 0.54 THz (Supplementary Sections X and XI). As shown in Fig. 4a, because of the redshift to the frequency ν red (0.53 THz), the SH and TH peaks are almost centered at 2ν red (1.06 THz) and 3ν red (1.58 THz), respectively.

Mechanisms responsible for the decay of the observed Faraday signals
The faster decay before 25 ps (Fig. 4c) of the fundamental oscillation amplitude at higher excitation intensities is attributed to two mechanisms. One mechanism is the magnetization-dependent Gilbert damping (Supplementary Section VI), which represents the decay rate of spin precession to the equilibrium state proportional to the magnetization change. Owing to this damping, a larger magnetization change can lead to faster decay. The second mechanism is considered to be the interference between oscillation components with different frequencies. When the probe pulses are transmitted through the sample, they experience different degrees of magnetization change because of the inhomogeneous distribution of the magnetic field along the z-axis (Fig. S4 in Supplementary Section II). Additionally, we need to consider that the frequency of spin precession shifts depending on the magnitude of the magnetization change. The interference among the components with different frequencies leads to a beating wave (Supplementary Sections XII and XIII) and thus results in a faster decay of the Faraday signal at higher excitation intensities.
Furthermore, the observation of an L-change-induced asymmetric oscillation in Δη motivates an investigation of the contributions of both M and L to the harmonic oscillations. Figure 4e shows the spectra of theoretically predicted Δη(t) results obtained under different assumptions. The blue and red curves are the spectra of the calculated results in Fig. 1d that reproduce the data for E THz = 0.2 and 0.8 MV cm −1 , respectively. On the other hand, the black dashed curve is the spectrum for E THz = 0.8 MV cm −1 obtained by fixing ΔL x /L 0 to zero. The comparison between the red and black dashed curves indicates that the observed SH peak at 2ν is caused not only by the second-order harmonic oscillation of ΔM z /M 0 , but also by the fundamental oscillation of ΔL x /L 0 (Supplementary Section VIII). This is because ΔL x /L 0 oscillates at twice the frequency of ΔM z /M 0 and thus contributes to the even-order harmonic oscillations. Meanwhile, the odd higher-order harmonics originate only from the nonlinear response of ΔM z /M 0 . Thus, it can be considered that the observation of the third harmonic is important: it helps us to differentiate between the nonlinearities of the M and L contributions, because the third harmonic peak is mainly determined by the oscillation of M z (as shown in Fig. S10 of Supplementary Section VIII). The second and fourth harmonic peaks contain significant contributions from both M z and L x . Thus, to obtain clearer experimental evidence that shows the nonlinearity of spin precession, the observation of the third harmonic spin oscillation is important.

Relation between the fundamental, SH and TH amplitudes
In HoFeO 3 , the Dzyaloshinskii-Moriya interaction leads to a canting of the sublattice magnetizations. Thus, the potential becomes anharmonic and the system has a broken symmetry, which allows even-and odd-order harmonic oscillations to be generated. To confirm that the observed harmonics in Fig. 4 originate from the nonlinearity of the spin response (and not from other effects such as a nonlinear response of the metallic structure), let us examine the dependence of the SH and TH signals on the fundamental oscillation. As shown in Fig. S16 (Supplementary Section XIV), the faster decay verified in Fig. 4c causes a deviation of the field-strength dependence of harmonics from the power law of the electric (magnetic) field. Each peak-area value plotted in Fig. S16 was obtained from the Fourier transform of the Δη waveform extending over 47 ps, and thus reflects both the amplitude and the lifetime of the fundamental, SH and TH components. As the excitation becomes stronger, the resulting faster decay of each component suppresses the increase in the corresponding peak area more effectively. This provides an intuitive explanation of the observed fieldstrength dependence. On the other hand, the SH and TH amplitudes, that is, the harmonics with order n = 2 and 3, closely follow the nth power of the fundamental peak, as shown in Fig. 5. The observed dependences are well reproduced by our calculations. These results indicate that, while the strong excitation condition modifies the fieldstrength dependence of the harmonics, the observed harmonics truly originate from a large magnetization change. In addition, our observation of the SH and TH is consistent with the selection rule derived from microscopic considerations on the canted antiferromagnet (Supplementary Sections XV and XVI).
In conclusion, owing to the strong excitation, harmonic oscillations up to the third order were observed. In particular, by realizing large magnetization changes, we were able to observe an L-changeinduced asymmetric oscillation. We have explained that the fundamental oscillation of ΔL x /L 0 contributes to the observed SH peak in addition to the second-order harmonic oscillation of ΔM z /M 0 , but it does not contribute to the TH because ΔL x /L 0 oscillates at 2ν. In addition, the ability to induce a large change in L is considered very important for future applications of ultrafast spintronics, because the absolute value of L is much larger than that of M, and thus the dynamics of L determine the spin current injected by spin precession 10,11,32 . Beyond providing an understanding of the spin nonlinearity and potential applications of antiferromagnets, our achievements indicate that our method may serve as a tool to control the functional properties of solids, including ferrimagnets, multiferroelectric materials, and quantum spin systems 18,33,34 , by strong THz magnetic near fields. . e Spectra of theoretically predicted Δη(t) results under different assumptions. The black curve is the spectrum for E THz = 0.8 MV cm −1 obtained in the case that ΔL x /L 0 is fixed at a value of zero, while the red dashed and blue dotted curves are, respectively, the spectra for E THz = 0.8 and 0.2 MV cm −1 obtained without fixing ΔL x /L 0 to zero, respectively. Note that "a.u." is the abbreviation for "arbitrary units".

Sample properties and spin dynamics in a HoFeO 3 crystal
We used a HoFeO 3 single crystal with a c-cut surface in the Pbnm setting 35 . The crystallographic a-, b-, and c-axes are parallel to the x-, y-, and z-axes, respectively. The HoFeO 3 single crystal was grown by using the floating-zone method and polished to a thickness of 52 µm. A single magnetic domain was confirmed by a two-dimensional static Kerr ellipticity measurement at zero magnetic field. Before each experiment, we applied a DC magnetic field to the sample to saturate its magnetization along the z-axis. The hysteresis in the Kerr rotation plotted versus the magnetic field is shown in Supplementary Section I. We fabricated an array of gold microstructures with a thickness of 200 nm on the crystal surface by an electron-beam lithographic process.
At room temperature, the two magnetizations m i (i = 1, 2) of the iron sublattices in HoFeO 3 are mainly antiferromagnetically aligned along the x-axis (with a slight canting angle β 0 = 0.63°due to the Dzyaloshinskii-Moriya field) and exhibit a spontaneous magnetization M along the z-axis (see the two gray arrows and the green arrow below the sample in Fig. 1a) 36 . In the THz region, there are two antiferromagnetic resonance modes, the so-called quasi-antiferromagnetic (q-AF) mode and the quasi-ferromagnetic (q-F) mode. The magnetic near field B z in our setup causes a q-AF-mode motion; as illustrated by the gray dashed curves in Fig. 1a, the Zeeman torque exerted on the spins triggers a precessional motion of each sublattice magnetization m i about the equilibrium direction. This precession causes an oscillation of the macroscopic magnetization M = m 1 + m 2 in the z-direction, as shown by the green double-headed arrow. The resultant magnetization change in the z-direction, ΔM z , is detected as the Faraday ellipticity signal Δη.

Faraday ellipticity measurements
In this experiment, we used a THz electric field to induce an oscillating current in the spiral-shaped structure and generate a strong magnetic near field B z perpendicular to the HoFeO 3 crystal surface. The magnetization change induced by the THz magnetic field is recorded by time-resolved measurements of the Faraday ellipticity. An amplified Ti:sapphire laser (repetition rate 1 kHz, central wavelength 800 nm, pulse duration 80 fs, and 7 mJ/pulse) was used to generate intense THz pulses by optical rectification of the femtosecond laser pulses in a LiNbO 3 crystal used in the tilted-pump-pulse-front scheme. The THz pulses were focused on the gold microstructure by using an off-axis parabolic mirror with a focal length of 50 mm, resulting in a spot diameter of ≈300 µm (full width at half-maximum). As shown in the inset of Fig. 1b, the spectrum of the incident-filtered THz electric field had a peak at around 0.57 THz. For the time-resolved Faraday ellipticity measurements, we used a 50× objective lens to focus linearly polarized optical probe pulses in the plane of the microstructure (spot diameter 1.2 µm), which enabled us to measure local changes in the Faraday ellipticity angle η. Note that the HoFeO 3 crystal was sufficiently transparent for the used 800-nm probe pulses. The polarization directions of the THz light and the probe light were parallel and almost along the y-axis. A balanced detector combined with a quarter-wave plate and a Wollaston prism was used to measure the ellipticity angle of the reflected probe pulse. The detection geometry and the definition of η are shown in Fig. S1(Supplementary Section I). All experiments in this study were performed at room temperature.

Data availability
Source data are available for this paper. All other data that support the plots within this paper and other findings of this study are available from the corresponding author upon request. Source data are provided with this paper.

Code availability
The code used to simulate the magnetization dynamics and the Faraday ellipticity change is available from the corresponding author upon request.